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Abstract. Using Molecular Dynamics (MD) and Monte Carlo (MC) simulations 
intcrfacial properties of crystal-fluid interfaces are investigated for the hard 
sphere system and the one-component metallic system Ni (the latter modeled 
by a potential of the embedded atom type). Different local order parameters 
are considered to obtain order parameter profiles for systems where the crystal 
phase is in coexistence with the fluid phase, separated by interfaces with (100) 
orientation of the crystal. From these profiles, the mean-squared interfacial width 
is extracted as a function of system size. We rationalize the prediction of 
capillary wave theory that ui^ diverges logarithmically with the lateral size of 
the system. We show that one can estimate the intcrfacial stiffness 7 from 
the interfacial broadening, obtaining 7 O.bkgT/cr^ for hard spheres and 
7 « 0.18 J/m^ for Ni. 



1. Introduction 

A key issue towards the microscopic understanding of crystallization from the melt 
is the knowledge about the properties of the equilibrium crystal-melt interface. 
Although experimental techniques such as electron microscopy and X-ray scattering 
give insight into the structure of solid- liquid interfaces [I], at least for atomistic 
systems, interfacial properties such as interfacial tensions or kinetic growth coefficients 
are hardly accessible in experiments. In principle, the situation is different for colloidal 
systems where microscopy allows for the direct measurement of particle trajectories 
and thus, similar as in a computer simulation, any quantity of interest can be computed 
from the positions of the particles. Recently, several experimental studies [51[31[31[S] 
were devoted to the study of solid-liquid interfaces in colloidal suspensions using 
confocal microscopy. However, a direct measurement of the anisotropic interfacial 
tension for a solid-liquid interface has not been realized so far. Moreover, it is an open 
question to what extent typical colloidal systems such as hard spheres can serve as 
model systems for crystallization processes on the atomistic scale, as they occur, e.g., 
in metallic alloys. 

Interfacial properties are also central parameters in the continuum modeling of 
crystal growth; e.g. the widely used phase field method needs the anisotropic interfacial 
tension as input (for a recent review of the phase field approach see Ref. [B]). The 
fact that interfacial tensions are in general not known from experiments reduces the 
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predictive power of the phase field method. Recently, more microscopic approaches for 
the description of crystallization processes have been proposed. The phase field crystal 
(PFC) method [7] is a generalization of phase field modeling to the atomistic scale. 
As shown by van Teeffelen et al. [5], the PFC method can be derived from dynamic 
density functional theory (DDFT) [9l [10] , the latter providing an "ab initio approach" 
to dynamic crystallization and freezing phenomena. Thus, there is hope that both 
PFC and DDFT will lead to some progress towards a microscopic understanding of 
crystallization phenomena. We note also that in the framework of static density 
functional theory, thermodynamic properties of the crystal-melt interface such as 
the interfacial tension can be predicted, at least for simple model systems (e.g. hard 
spheres) pT] . 

For all the latter models, one has to keep in mind that they introduce various 
severe approximations: in particular, since they are mean-field theories, statistical 
fluctuations are neglected. Thus, capillary wave excitations [12] are not taken into 
account that strongly affect the interfacial properties (e.g. a broadening of the mean- 
squared interfacial width). 

Beyond mean-field theories, particle-based computer simulations are an 
appropriate tool to study solid-liquid interfaces at a microscopic level. In principle, 
molecular dynamics (MD) as well as Monte Carlo (MC) simulations provide a 
numerically exact treatment of the statistical mechanics, only based on a model 
potential that describes the interactions between the particles. However, when 
one examines the details of the solid-liquid interface at phase coexistence via such 
simulation techniques, one must be aware of various problems: a slight deviation 
from the correct coexistence condition, details of the averaging procedure, insufhcient 
sampling of statistical fluctuations, or an inappropriate choice of the order parameter 
may cause more or less drastic deviations from the correct result. These problems have 
hampered progress in the area. On the one hand, new and presumably rather accurate 
methods for the analysis of solid-liquid interfaces have been presented for hard spheres 
[131 fT4] , Lennard- Jones systems [15] , models of metallic systems [Ml [171 Ull [12] ; and 
silicon [5D]. On the other hand, in none of these studies it has been systematically 
analyzed how finite size effects influence the properties of solid-liquid interfaces. In 
the present work, we do the first steps to fill this gap and we demonstrate that, in the 
framework of capillary wave theory, the interfacial stiffness can be estimated from an 
analysis of finite size effects. 

A comparative study is performed of the structure of solid-liquid interfaces of hard 
spheres and a model of Ni. In this manner, we shed light on differences between both 
kinds of systems and thus, we contribute to the question whether a typical colloidal 
system (hard spheres) can serve as a model for a typical metallic system (Ni) with 
respect to the interfacial properties at coexistence. Furthermore, by using both MC 
and MD simulations we rationalize that both methods can be applied to the interfacial 
analysis. 

The outline of the paper is as follows. In the next section, we introduce the 
models used for the simulations and summarize the simulation methodology. In Sec. 3 
we give a brief overview over the results of capillary wave theory that we use to 
determine the interfacial stiffness. Then, in Sec. 4 we present the results for the 
fluid structure, the local order parameter profiles and the system size dependence of 
the mean-squared interfacial width from which we estimate the interfacial stiffness. 
Finally, we summarize the results. 
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2. Models and simulation techniques 

In this section, we introduce briefly the main details of the MC simulations for the 
hard sphere system and the MD simulations for Ni. Moreover, for both cases we 
describe how we have prepared configurations in a slab geometry where the crystal 
phase in the middle of an elongated (rectangular) simulation box is at coexistence 
with the fluid phase, separated by two interfaces (parallel to the xy-plane) and using 
periodic boundary conditions in all three spatial directions. 

2.1. Monte Carlo simulation of hard spheres 

The hard sphere system is defined via the potential 



with r the distance between two particles and a the diameter of a particle. Throughout 
the paper all length scales are measured in units of a for the hard sphere system. 

Since for any allowed hard sphere configuration, the total potential energy is 
zero, temperature T is only a scaling factor and the thermodynamic properties are 
fully controlled by the packing density 77 = (or, when the total volume V 

of the system is a fluctuating variable, by the pressure P). As a consequence, the 
phase behavior of hard sphere systems is completely driven by entropy |21] . As first 
claimed by Kirkwood [H], hard spheres exhibit a fluid-to-solid transition. However, 
Kirkwood's prediction was based on misleading arguments and so it was a surprising 
discovery when the freezing of hard spheres was first observed in early computer 
simulations [23l[2l]. Since the phase behavior of hard spheres depends only on packing 
density 77, the phase diagram is particularly simple. Due to the absence of attractive 
interactions, the hard sphere model ([T|) does not exhibit a liquid-vapor transition. 
Fluid and fee crystal coexist between the freezing point ?/f = 0.494 and the melting 
point ?/in = 0.545, while the pure crystal is the stable phase for 77 > 77m [251127]. 

The MC simulations were carried out in the isothermal-isobaric {NPT) and 
(NPzT) ensemble, in which the pressure P, the temperature T and the number of 
particles TV are constant (in the NP^T ensemble, P^ is the diagonal component of the 
pressure tensor perpendicular to the xy plane) . A MC code was developed applying the 
standard Metropolis algorithm [26l [27l [28] . The trial moves are particle displacements, 
where each particle was attempted to be displaced once per MC cycle, and system's 
volume rescaling that was executed once per MC cycle. The maximum displacement 
is chosen to keep the acceptance rate at 30% for the particles and 10% for the volume. 

In test runs, we have reproduced the solid and liquid branches of the equation 
of state and have found full agreement with analytical results [29l [30] . Additionally, 
the radial distribution function and the static structure factor for the bulk liquid has 
been compared to the Percus-Yevick approximation [31]. The coexistence pressure 
was found by the interface velocity analysis similar to Ref. [Hj. At coexistence the 
total volume of the melt-crystal system is constant, since the melt is in equilibrium 
with the crystal, therefore the interface velocity is vi = 0. Hence, a series of runs 
were performed in a wide range of pressures. At each simulation, the interface 
velocity vi was estimated from the slope dV/dt of the temporal change of the 
system volume V{t) in the stationary growth interval. Following this procedure, 
first, we have estimated the rough location of the coexistence pressure and then 
increasing the resolution of the pressure interval we improved the initial accuracy. 




(1) 
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The final result yielded = 11.55 ± 0.04 kgT/a^, very close to the literature value 
Pc = 11.567 keT/a^^. Similar values {pc = 11.54,11.5,11.53,11.55 keT/a^) have 
been reported in Ref. [131 1311 1331 131| accordingly. The average densities were found 
to be ps = 1.04 for the solid and p\ = 0.9385 for the liquid, close to the results 
of Hoover and Ree [25] {ps ~ 1.04 for solid and p\ ~ 0.939). The bulk interplanar 
distance between lattice planes in the crystalline phase is ~ 0.784 <t. Further details 
are discussed elsewhere |35j . 

To generate crystal-melt interfaces at coexistence, we first prepared solid-liquid 
"sandwiches" where the (100) direction of the fee phase is oriented perpendicular to 
the z axis. In this study, we considered the system sizes of side length L = na lattice 
spacings with n = 5, 6, 7, 8, 10 (where a = 1.567ct at coexistence). The total number 
of particles is iV = 2500, 4320, 6860, 10240, 14580, and 20000, respectively. We 
have verified that the elongation = 5L along the z direction is sufficient to avoid 
interactions between the interfaces due to periodic boundary conditions. For smaller 
elongations in z direction we see significant finite size effects in density profiles due to 
the interaction of the two interfaces via periodic boundary conditions. 

As a first step, a solid slab of size {LxLx 2>L) and a liquid box of size {LxLx 2L) 
were equilibrated separately for 10^ MC cycles in the NPT ensemble at Pc- Here, L 
corresponds to the bulk solid density ps at coexistence. For the liquid, the same values 
of L were used in x and y direction. Then, the correct bulk density of the liquid was 
achieved by using simulations in the NP^T ensemble. Next, the two parts were placed 
together in a simulation box of size x Ly x with ~ Ly ^ L. 

The most delicate point of the preparation is how to match the solid and fluid 
parts. The initial configuration might be refused if the melt and the crystal are 
too close, otherwise a large gap between them would cause lower fluid density and 
artifacts with respect to interface properties. Therefore, new fluid-solid configurations 
were relaxed in NP^T simulations until the coexistence densities were recovered in 
the bulk regions. The first 10^ MC cycles the lateral sizes and the positions of the 
solid particles were fixed in order to avoid internal stresses in the solid. Furthermore, 
for each system size we have performed isochoric runs of 2 • 10^ MC cycles, initially 
rescaling the length of the box to the average value {L^), as computed in the last 
5 • lO'* MC cycles of previous NP^T runs. Next, we have selected 50 independent 
configurations for sizes n = 5, 6, 7, 8 and 10 independent configurations for n = 10 
to compute the equilibrium properties of the interfaces over 5 • lO'' MC cycles in the 
isochoric ensemble. The statistics was collected every 20 MC cycles. The length of 
the final runs was chosen to prevent the diffusive motion of the interface, however 
in several configurations an additional interface broadening was detected due to the 
displacement of the solid-liquid interface. Those runs were replaced by additional 
runs. 

2.2. Molecular Dynamics simulation of Ni 

For Ni, a similar methodology as for the hard sphere system has been employed. But 
MD simulations instead of MC simulations were performed. As a potential to model 
the interactions between the particles in Ni we used a potential of the embedded atom 
type, as proposed by Foiles [36]. We show elsewhere [37] that this model potential 
reproduces very well various thermodynamic and transport properties, as obtained 
from experiments of liquid Ni. 

As before for the hard sphere system, an inhomogeneous solid-liquid system is 
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simulated in a sandwich geometry with L x L x being the size of the simulation 
box. Again, = 5L is chosen. In the following, L will be also expressed in terms 
of the number of lattice planes of the fee crystal, n, as L = na (note that the lattice 
constant is given by a = 3.58 A at the melting temperature T^). 

Newton's equations of motion are integrated with the velocity form of the Verlet 
algorithm [38], using a time step of Ifs. The melting temperature of the model is 
again estimated from an interface velocity analysis, considering an inhomogeneous 
solid-liquid system in the NPT ensemble. The melting temperature is obtained from 
the linear fit of the interface velocity versus temperature up to an undercooling of 
about 40 K. From these simulations [37], we have found the melting temperature 
Tm = 1748 K, which is in good agreement with the experimental value, = 1726 K 
[39] . Simulations of different system sizes indicate that finite size effects are weak [37] , 
as far as the determination of is concerned. For the Ni model used in this work, the 
melting temperature of the smallest system with N = 2500 particles was 0.5% higher 
than the estimated melting temperature in the thermodynamic limit, Tm = 1748 K. 

To prepare an inhomogeneous system with two crystal-liquid interfaces at Tm, 
the following steps are involved: First, atoms, disposed on a fee lattice, are relaxed 
in a NPT simulation for about 30 ps at Tm and zero pressure. Temperature was kept 
constant by coupling the system to a stochastic heat bath, i.e. by reassigning every 200 
steps new velocities to each particle according to a Maxwell-Boltzmann distribution. 
To keep the pressure constant, an algorithm proposed by Andersen was used, setting 
the mass of the piston to lOOeVps^/A^ [40]. In the next step, a liquid and a crystal 
region are defined in the system such that the crystal region in the middle of the 
elongated simulation box occupies a volume of L^. Atoms in the crystal region remain 
at fixed positions while the rest of the system is heated up to 2400 K which is well above 
the melting point. At this step, only volume changes with respect to the expansion 
and compression of the box in the direction perpendicular to the liquid-solid interface 
are applied, i.e. the simulation is done in the NPzT ensemble. In order to completely 
melt the liquid regions, the simulation runs over about 100 ps. Then, the temperature 
of the melt is set back to the initial temperature at which the crystal was prepared. 
A run over 50 ps in the NPzT ensemble is added where all the particles are allowed 
to move. From the last 10 ps of this run the average length of the box in z direction, 
(Lz), is determined. After the length of the simulation box in z direction is rescaled to 
(Lz), the simulation continues with a run over 20 ps in the NVT ensemble. From the 
last 10 ps of this run, the average total energy of the system is computed. The system 
is set to this energy by rescaling the velocities of the particles appropriately. Finally, 
microcanonical production runs over 1 ns are done from which the information about 
the interfacial properties are obtained. We considered systems of lateral size L — na 
with n = 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, and 15. The total number of particles in these 
systems is iV = 2500, 4320, 6860, 10240, 14580, 20000, 26620, 34560, 43940, 54880, 
and 67500, respectively. For each system size, 5 independent runs were performed. 

3. Capillary fluctuations 

In this work, we consider systems where a crystal phase is at coexistence with a 
fluid phase and the two phases are separated from each other by interfaces (two 
interfaces are formed due to periodic boundary conditions). As a matter of fact, the 
presence of the interfaces breaks the translational invariance, which is a continuous 
symmetry property of the underlying Hamiltonian. The latter leads to the occurrence 
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of Goldstonc excitations: long-wavelength transverse excitations, known as capillary 
waves, appear that are thermally driven undulations of the interface. At infinite 
wavelength (i.e. in the limit of wavenumber g — > 0), they describe an overall 
translational motion of the interface with zero energy cost. In the framework of mean- 
field approaches such as phase-field modeling, capillary wave excitations are neglected, 
although they strongly affect the interfacial properties. 

This can be seen in the framework of capillary wave theory (CWT) [12]. This 
theory describes the free energy cost AF of long-wavelength undulations of an 
interface. For three-dimensional systems, CWT predicts a logarithmic divergence 
of the mean squared width, w^, of the interface with the lateral system size L. As 
we see below, in a computer simulation, this result of CWT can be used as a method 
to determine the interfacial tension by measuring the mean-squared width w'^ for 
different system sizes. Whereas this method has been successfully applied to the 
Ising model [HI \^ [33] , polymer mixtures [JH [351 HS] , and liquid- vapor interfaces 
in the Asakura-Oosawa model for colloid-polymer mixtures [171I3H], it has not been 
used as a method to estimate the interfacial tension of crystal-fluid interfaces. In 
the latter case, many of the recent simulation studies on hard spheres (e.g. [13j). and 
metallic systems (e.g. [17]) have used an analysis of the capillary wave spectrum to 
determine the interfacial tension. However, here, we demonstrate that also in the case 
of solid-liquid interfaces the interfacial stiffness can be computed by measuring as 
a function of InL. Both for polymer mixtures [44[ 145] and the liquid- vapor transition 
of the Asakura-Oosawa model [48] , it has been shown that the analysis of the capillary 
wave spectrum and the finite-size analysis of the interfacial broadening yield values 
for the interfacial tension that are in agreement. 

As before, we consider atomically rough crystal-fluid interfaces. We parametrize 
the local fluctuations of the interface by a function h{x, y) that denotes the local 
deviation of the interface position zo(x,?/) from the mean value h(x,y) = ZQ{x,y) — 
{zo{x,y)) (here, z is the Cartesian component perpendicular to the interface, x, y 
the ones parallel to it). h{p) = h{x,y) can be expressed in Fourier coordinates, 

= ^^h{q) e'>cp{iq • p) (with the wavevector q = {qx,qy)). Then, the total free 
energy of the interface can be expressed in reciprocal space as 

^^^^E^'Kf- (2) 

Here, is the area of the flat interface. 7 is the interfacial stiffness, defined by 
7 = 7-1- d'^j/dO'^ with 7 the interfacial tension and 9 the angle between the interface 
normal and the (100) direction. The interfacial stiffness takes into account the 
anisotropy of the interfacial tension in case of a crystal-fluid interface; of course, in 
case of a liquid- vapor interface 7 would be replaced by 7 in Eq. ([2]) . 

Since the different q modes in Eq. ^ are decoupled, it follows from the 
equipartition theorem that each mode carries an energy of ksT and thus one obtains: 

L'^'jq'^ 

This expression can be used to determine 7 by measuring the slope of the straight line 
that fits l/(|/i^p) plotted as a function of q^. In fact, this has been done in recent 
simulation studies of hard spheres [55] and metallic systems [T7]. However, in the latter 
works, a geometry with Ly << were chosen such that only the fluctuations of a 
quasi one-dimensional "ribbonlike" interface are considered. Although this simplifies 
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the analysis it also alters the nature of the capillary fluctuations (see below). In this 
work, we consider therefore only a geometry with L ^ L^. ~ Ly. 

Equation ^ can be used to determine the mean-squared interfacial width w^, 
given by 



{\hp\') = E(IM') = (1^/ dq{\h,f) , (4) 



which yields 



2 ksT f'^/'dq ksT 
w'cw = ir^ / — = ^ In . (5) 
27r7 J2^/L 1 27r7 

In Eq. ([5|) , £ is a cut-off length that is introduced in accordance with the assumption 
that only modes with a wavelength larger than the typical width of the interface 
are taken into account. Note that with the aforementioned quasi one-dimensional 
ribbonlike interface the mean-squared width would increase linearly with L. 

In mean-field theory the interface between coexisting phases is assumed to be flat 
and the interfacial profile (/'(z) is described by a hyperbolic tangent, 

Mz) = A + B tanh ( (6) 
\ J 

where A and B arc parameters related to the bulk values of the densities or order 
parameters, zq and wo are the position of the interface and its width, respectively. 

Now, the idea is to combine the mean-field result with that of CWT by considering 
wq as the width of an intrinsic profile that is superposed by fiuctuations described by 
w^^. The total width of the profile is then obtained from a convolution approximation 

2 2,^2 2 , 1 r 

w =w^ + -w^^ = Wo + ^^\nL- ^^Xni . (7) 

When using this equation as a fit formula to analyse profiles as measured in the 
simulation, it is impossible to disentangle the intrinsic width contribution iuq from the 
"cut-off" contribution — In^. This issue has been discussed in detail for the case 
of polymer mixtures [H SSI Si [501 [H] . 

The main issue in the following is to investigate whether Eq. ([7]) can be used 
to measure the interfacial stiffness in a computer simulation. To this end, fits with 
Eq. ^ to order parameter profiles are used to obtain an effective interface width for 
different lateral system sizes L. This will be described in the next section. 

4. Results 

^.1. Static structure factor: Hard sphere fluid vs. Ni melt 

In order to study the structural differences between the bulk hard sphere fluid and the 
bulk Ni melt at coexistence. Fig. [1] displays the static structure factor [31] for both 
systems, 

S{q) = ]^l^Y.GM^q-rk]^^ (8) 

with r/c the position of particle k and q the three-dimensional wavevector (different 
from the two-dimensional wavevectors that we have considered in the previous section). 
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Figure 1. Static structure factor S{q) of the iiard sphere system and the Ni melt 
at coexistence, i.e. at the volume packing fraction rj = 0.492 for the hard spheres 
and at the temperature T = 1748 K for Ni. 




Figure 2. Density profiles for a) hard spheres and b) Ni. In each case, the profiles 
are shown for two different system sizes. 



To provide a better comparison between the structure factors for the two systems, 
we have multiphed the wave-number q in Fig. [l] by ct' ~ a for the hard sphere 
system and a' = 2.24 A for Ni (which is similar to the nearest neighbor distance, 
rNiNi =2.42A [37]). 

Mainly two differences between the two systems can be inferred from Fig. [T] 
Towards q —>■ the structure factor for Ni has a much lower amplitude, indicating 
that, at coexistence the Ni melt has a lower compressibility than the hard sphere 
system. Moreover, the amplitude of the first peak in S{q) for Ni is slightly larger. But 
the overall shape of S{q) is surprisingly similar for both systems which shows that the 
hard sphere system's structure is very similar to that of Ni. 
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Figure 3. Probability distributions of the order parameter in Ni (broken lines) 
and the hard sphere system (solid lines) for the liquid and the fee phase, as 
indicated; a) qsqe, b) ^. Different time intervals for obtaining the average particle 
positions are considered for Ni; thereby, the dashed lines correspond to an interval 
of 0.1 ps and the dotted lines to 1 ps (the latter case is referred to as Ni*). 



4-2. The structure of the solid-fluid interfaces 

A simple quantity to characterize the structure of the interface is provided by the 
density profile p{z) across the interface. To compute p{z), the system is divided into 
slices of thickness Az and then, in each slice, one counts the number of particles and 
divides by the volume of the slice L^Az. The displacement of the lattice planes along 
the z-axis during the time evolution was corrected for each configuration. As for the 
order parameter profiles shown below, we have averaged the density profiles over the 
two interfaces in each system that are present due to periodic boundary conditions. 

Figure [2] shows density profiles for the hard sphere system and Ni, in each case 
for two system sizes. The shape of the profiles is typical for inhomogeneous systems 
with crystal-fluid interfaces. In fact, very similar density profiles have been found for 
various systems with solid-liquid interfaces [13l [151 Ull HH HO] • Whereas one observes 
huge oscillations in the crystalline region due to the presence of crystalline layers, in 
the fluid region p(z) is constant. In between, i.e. in the interface region, the amplitude 
of the peaks decreases. Obviously, the size effects are very small for both considered 
systems. However, for the big systems the height of the peaks is slightly larger in the 
interface region. We note that although the density proflles have a very different shape 
in the solid and the liquid regions, both for the HS system and for Ni the differences 
in the average densities of the solid and the liquid phase are rather small. In the HS 
case, the solid density is about 10% higher (ps = 1.04, pi — 0.9385) and for Ni, it is 
about 5% higher {ps = 8.357g/cm'^, p\ = 7.928 g/cm"^). Due to this small difference 
between pi and ps the density is not a good order parameter for the investigation of 
interfacial properties such as the interfacial width. 

Another possibility to characterize the structure of interfaces is provided by 
proflles of local order parameters. Steinhardt et al. [Slj have proposed rotational- 
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invariant order parameters in terms of expansions into spherical harmonics Yim, 



with 

1 ^' 

Qlm{i) = ^ E '^(^^■)) , (10) 

* j = l 

where rij is the distance vector between a pair of neighboring particles i and j, is 
the number of neighbors within a given cut-off radius, and 0{rij) and (t>{rij) are the 
polar bond angles with respect to an arbitrary reference frame. 

Similar local order parameters have been introduced by ten Wolde et al. [53] . 
defined by 

1 

qmii) ^^^liiii) ■ mU) ■ (n) 

The internal product in this sum is given by 

qi(i)-qi(j)= E ^im{i)mm{jY (12) 



f ■\ Qlm{i) /To\ 
llmW = — ■ (13) 



with 



In the following, we use the parameter qeqe that is defined by Eqs. (fTT|) - (fT3|) . setting 
/ = 6. 

Another local order parameter used in this work was introduced by Morris |54j 
11^-^' 2 
1 ' j=i k=i 

where the wavevectors q^ are chosen such that in a perfect crystal 

I exp(igfc • fij)\ = 1 . (15) 

Again, r*y is the distance vector between neighboring particles. With respect 
to the basis vectors of the fee lattice with lattice constant a, oi = a/2(l,l,0), 
a2 = a/2(0, 1,1) and 03 = a/2(l,0, 1), an appropriate choice of wavevectors is 
61 = 27r/a(-l,l,-l), 62 = 27r/a(l,-l,l) and 63 = 27r/a(l, 1, -1). An additional 
average of ^'(i) over a particle with index i and its neighboring particles yields 

(16) 

The parameter 4* together with q^qe is used in the following to distinguish solid 
particles from fluid particles and to identify the interfacial regions. To select the 
nearest neighbors we introduced the cut-off radii that correspond to the first minimum 
of the radial distribution function of the bulk liquid phase at coexistence. 
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Figure 4. Profile for the order parameter 9696, a) for hard spheres and b) for 
Ni. The inserts provide a magnification of the interface region. 




Figure 5. Profile for the order parameter 'I'r, a) for hard spheres and b) for Ni. 
The inserts provide a magnification of the interface region. 



Figure [3] displays local order parameter distributions for the pure liquid and the 
pure solid phases. The distributions indicate that the considered order parameters 
are well-suited to distinguish liquid from solid particles. For the calculation of the 
local order parameters, we used time-averaged particle positions. To this end. for the 
hard sphere system, positions were averaged over 50 MC cycles. For Ni, the phonon's 
degrees of freedom lead to a significant shift of the order parameter distributions. 
Using a time interval of 0.1 ps for the averaging of the particle positions in Ni, compared 
to an ideal fee crystal the order parameter distributions for the crystal phase are 
broader and tend to shift to smaller values of the order parameter (see Fig. [S]). For 
a time average over 1 ps the order parameter distributions are very similar to those 
for the hard spheres. However, since in Ni at Tm a time scale of 1 ps is already close 
to the time scale of particle diffusion in the melt, we have used a time averaging over 
0.1 ps for the analysis of the intcrfacial properties that are presented in the following. 

In Figs. [4] and [5l profiles of the local order parameters are shown, i.e. the sum 
of the order parameter of the particles contained within a slice transversal to the 
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Figure 7. Mean-squared width as a function of InL, a) for hard spheres and 
b) for Ni. The values for 7 are obtained from the fits (solid lines). 



solid-liquid interface divided by the volume of the slice Ay ~ L^/S.z. The profiles are 
calculated with a resolution (bin size) of 0.02(7 for the hard spheres and 0.1 A for Ni. 
The order parameter profiles show similar features as the density profiles. However, 
finite-size effects seem to be revealed in a more pronounced manner. Clearly, the 
height of the peaks in the interface region slightly increases with increasing system 
size. 

To compute coarse-grained profiles we identify the minima in the fine-grid profiles 
of Figs. [4] and [5l These minima define the borders of nonuniform bins that match the 
crystalline layers. Then, we compute the average value of the order parameter in 
each of the latter bins. Examples for the resulting coarse-grained order parameter 
profiles for the case of q^q^ are displayed in Fig. [S) Here, the solid lines are fits with 
a hyperbolic tangent function, 0(2;) — A — i?tanh[(z — zq)/w\ (where A and B are 
parameters related to the bulk values of the order parameter, and zq and w are the 
interface position and its effective width, respectively). Both for the hard spheres and 
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Ni. the latter fits indicate that the width w is larger for the big systems, as expected 
from CWT. 

4-3. Estimate of the interfacial stiffness 

In Fig. [21 the mean-squared width w^, as obtained from the fits to the coarse-grained 
profiles for q^qe and ^, is plotted as a function of InL. The plot confirms the 
logarithmic increase of w'^ with the lateral system size, as predicted by CWT. 

From the fits with Eq. ([7]). we estimate for the (100) orientation 7 ~ 0.50 ± 
0.05 /a^ for the hard-sphere system and 7 ~ 0.18±0.01J/m'^ for Ni. These values 
roughly agree with previous estimates, obtained by other methods. For hard spheres, 
Davidchack and Laird [14] found 7 = 0.57 ksT/a^ using a thermodynamic integration 
approach, while Mu et al. [55] obtained 7 c± 0.62 kBT/a^ from the analysis of the 
capillary wave spectrum. However, Davidchack et al. |56| later criticized their result 
as being biased and rather suggested 7 — 0.56 ± 0.02 kBT/a^ when averaged over all 
interface orientations. For the (100) orientation, they suggest 7 ~ 0.44±0.03 ksT/a^ . 
However, their actual data reveal huge fluctuations, and the judgement of the actual 
accuracy may need reanalysis. For Ni, Hoyt et al. [17| determined the interfacial 
stiffness for the (100) orientation (and other orientations of the Ni fee phase). Using 
a different EAM model and the analysis of the capillary wave spectrum to measure 7, 
they obtained 7 ~ 0.23 J/m^, which is slightly larger than our result. 

5. Summary and Outlook 

In this paper, we have presented a comparative study of melt-crystal interfaces for hard 
spheres and an embedded atom model for nickel. These rather diverse systems have 
been studied in analogous geometries, namely L x L x rectangular simulation boxes 
with periodic boundary conditions, at conditions where a crystalline slab, separated 
by two L X L interfaces oriented perpendicular to the z-direction, coexists with the 
fluid phase. The motivation for this study was to provide a better understanding of 
the information that one can extract from the simulation study of such interfaces, 
paying particular attention to finite size effects, and to limitations of the accuracy 
which are inherently due to the simulation setup. In fact, due to the translational 
invariance of the simulation geometry as a whole, the center of mass of the crystalline 
part is not fixed in space, but may fluctuate and diffuse along the z-axis. In addition, 
at phase coexistence in a finite box, fluctuations may occur when the size of the crystal 
(volume fraction of the box that is crystallized) changes. As a consequence, in each 
configuration that is analyzed one must locate the precise position along the z-axis 
where the lattice planes are (this was done via a fine-grid coarse graining) and then the 
time averages are found in such a way that the positions of lattice planes and interface 
centers coincide. It is clear that this is a delicate procedure and hence there is the 
need to watch out for possible systematic errors, which are not necessarily equally 
important for different kinds of systems, and for different simulation methods (e.g. 
MC and MD). In view of these caveats, it is gratifying to state that with the methods 
described in this paper, these problems seem reasonably well under control. Indeed, 
we find that the main difficulty in the interpretation of our results for the interfacial 
profiles and their width is the broadening by the capillary waves. This phenomenon, 
though well-known for vapor-liquid interfaces, has found comparatively little attention 
for the melt-crystal interface. Our results imply that the capillary wave broadening 
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(for rough, non-facetted crystal surfaces) is present and important: while it makes 
a naive direct comparison with DFT calculations of interfacial profiles obsolete, it 
yields a relatively straightforward method for extracting the interfacial stiffness and 
the accuracy of this method seems to be competitive to other approaches. 

As a next step wc plan a detailed comparison with the alternative method where 
the Fourier spectrum of interfacial fluctuations is analyzed. For vapor-liquid type 
systems, such comparisons can be found in the literature, but for the melt-crystal 
interface a comprehensive comparative assessment of different methods still is lacking. 
Of course, for applications in crystal nucleation and growth phenomena very accurate 
estimates for the interfacial stiffness are indispensable. 
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